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Predicting and controlling the dynamics of infectious diseases 


Robin J. Evans 1 and Musa Mammadov 1 


Abstract —This paper introduces a new optimal control 
model to describe and control the dynamics of infectious 
diseases. In the present model, the average time of isolation 
(i.e. hospitalization) of infectious population is the main time- 
dependent parameter that defines the spread of infection. All 
the preventive measures aim to decrease the average time of 
isolation under given constraints. 

I. INTRODUCTION 

The outbreak of Ebola Vims Disease (EVD) in West 
Africa 2014 revealed many challenges in predicting and 
controlling the spread of infectious diseases. These chal¬ 
lenges are partly related to the mathematical modeling of 
the dynamics of the epidemic. Providing accurate predictions 
appeared to be extremely difficult. 

To address these challenges, several new models have been 
suggested each providing quite different results, for example 
[1], [2], [3], [4], [13], [14], [15], We also note that different 
aspects of possible control have been intensively studied in 
the literature including distribution strategies for vaccination 
and antibiotic programs [12] as well as travel restrictions [8], 

In this paper we concentrate on the development of models 
that are well suited to the control of an outbreak. The most 
commonly studied models in this area deal with temporal 
networks [10], [11], [20] where several different models have 
been suggested. 

In [9] the 1995 Ebola outbreak in Congo is considered us¬ 
ing an SEIR model whereby control intervention, performed 
at time t*, is described by the transmission coefficient /3(f) 
defined by /3(f) = /3 if f < f* and /3(f) = /3 exp(— q(t — f*)) 
for f > f*, where /3 is the initial transmission rate that would 
remain stable without the intervention. 

Our paper models and studies an alternative control mech¬ 
anism to decrease in the transmission of infection based on 
the model developed in [7] where the transmission of infec¬ 
tion depends mainly on two parameters: the transmission rate 
/3 and the average time of isolation r. In contrast to [9] we 
assume that the transmission rate /3 does not change over 
the whole period under consideration. The main dynamic 
parameter in our model is r and all the intervention measures 
are directed at decreasing r and consequently reducing the 
spread of infection. 

The average time to hospitalization can be used for the 
average time of isolation r, although the isolation of infec¬ 
tious population is not exactly the same as “hospitalization”. 
Clearly, as the disease progresses hospitals become short on 
beds (as well as staff and supplies) required to isolate and 
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treat all newly infected individuals. As a result, the number 
of infected population can grow exponentially. This was the 
case for the Ebola virus epidemic in West Africa (Guinea, 
Sierra Leone and Liberia) where the spread of infection 
was highly dangerous during June-November 2014 when the 
capacity for treating Ebola patients was insufficient. It was 
reported [5] that dining this period “... many clinics and 
hospitals in all three of the countries worst hit by Ebola 
have effectively been shut down”. 

After rapidly building new infrastructure and increasing 
the capacity of beds the outbreak slowed down significantly. 
Starting from January 2015, the epidemic has moved to 
the ending phase that involves ensuring “capacity for case 
finding, case management, safe burials and community en¬ 
gagement” ([17] - WHO, Ebola Situation Report, 28 Jan 
2015). Note that in [6] the hospitalization rate was the 
parameter showing the greatest change. 

Addressing these issues, this paper suggests new mathe¬ 
matical models that can be used to increase the efficiency 
of available resources. The main goal here is to keep 
the model as simple as possible and, at the same time, 
to have measurable control variables. Note that there are 
many useful control measures that have been intensively 
studied by introducing more “detailed” mathematical models, 
however such models have less predictive capabilities (due 
to overfitting). Prediction is crucial when considering future 
planning periods. 

The main component in the suggested model is the optimal 
distribution of bed capabilities across countries/regions. This 
is a very important and difficult problem that requires an 
accurate prediction of the dynamics of infected population 
in each region. For example, evaluating the situation of 
the Ebola outbreak, WHO’s Ebola Situation Report on 14 
Jan 2015 [17] notes that “Each of the intense-transmission 
countries has sufficient capacity to isolate and treat patients, 
with more than 2 treatment beds per reported confirmed and 
probable case. However, the uneven geographical distribution 
of beds and cases, and the under-reporting of cases, means 
that not all EVD cases are isolated in several areas.” 

II. MODEL 

In [7] a new model is introduced to study the dynamics 
of epidemics by considering the average time for isolation 
(denoted by r) of infectious population as a time-dependent 
parameter. This model is derived from the well studied SIR 
(Susceptible-Infectious-Recovery) model (e.g. [16]) and is 
similar to models based on transmission rates from infectious 
population at different generations (e.g. [13]). 



The use of time-dependent parameter r enables the analy¬ 
sis of future scenarios by considering possible changes in r. 
In this paper we extend this approach by developing practical 
and efficient optimal control models. 

We denote by x(t) the number of infected cases at i £ 
{1, 2, ■ ■ •, T} (in days). Assuming that the natural death rate 
of population (p) is zero, the equation for x(t) is as follows 
(see [7] for more details) 

r—1 

x(t + 1) = j3 ^^(1 — au>(i)) x(t — d — i). (1) 

i=0 

Here a is the death rate due to disease; d is the average 
latent period (in days) for infected individuals to become 
infectious; r is the average infectiousness period (in days); 
it is the average time required for isolation (time to hospi¬ 
talization); and /3 is the transmission rate. Moreover, u is a 
gamma (cumulative) distribution function (with p.d.f - oj p ) 
for deaths due to disease. The fraction (1 — au(i)) in this 
case represents the proportion of remaining infected cases 
x{t — d — i) after d + i days. 

The sum 

T— 1 

I a (t) = ^(1 “ auj (i)) x(t-d - i) 
i =o 

defines the number of ’’active” infectious population at time 
f; it represents the number of infectious population that 
are not yet isolated and therefore it is the only source of 
secondary infections, (for the sake of simplicity we do not 
consider infections in hospitals and death ceremonies). Then 
by setting x(t) = (3I a {t) we obtain model (1) in [7] where 
// (the natural death rate) in our case is zero. 

The basic reproduction number R is calculated by consid¬ 
ering the stationary states in ([T}: 

T — 1 

R = f}[T-a'y' j u(i)\. (2) 

2—0 

There are three main parameters in (|T]> - a,f3 and r. The 
results obtained in [7] show that this model provides quite 
good approximation to the total number infected cases and 
deaths during the current Ebola epidemic if r is a piecewise 
constant function (in fact, constant over consequent subin¬ 
tervals with durations 2-3 months) and the parameters a and 
P are constant over the whole period. 

These results help us to predict the dynamics of an infected 
population at future time intervals by keeping the values of 
a and /3 unchanged (estimated from the previous period) and 
considering different possible changes in r. In this case the 
major strategy of preventive intervention is the achievement 
of some decrease in r that according to (|2]) is equivalent to 
decreasing the effective reproduction number. 

This approach is implemented below by introducing 
an optimal control models where the average time to 
hospitalization r is the key variable. According to the 
results of data fitting mentioned above, it is sufficient to let 


r be constant on quite long time intervals (months). 

Control t. Therefore, we define r(t) as a control 
variable by assuming that it is piece-wise constant with 
integer values (days). For the sake of simplicity let 

r(t) = TiGU, Vf e (T jt T j+1 ], J = 1,2,--- ,p. 

It is reasonable to assume that U = {T rn i n , T rn i n + 
1)'''! T max}'i where T m i n > 1 is the minimal number of 
days required to isolate infectious population. 

Trajectory x. Given control r we define trajectory 
x = x(t) as follows 

r(i)-l 

x(t + 1) = /? (1 — x(t — d — i). (3) 

»=o 

In this formula the sum X^I=d _1 (l — aui(i))x(t — d — i) 

represents the number of infectious individuals that are not 
yet isolated. Considering the average length of hospital stay 
(in days), the number of hospitalized cases at t can be 
calculated as 

a 

h(t : r, x) = ^ (1 — aui(i)) x(t — d — i) (4) 

i=r(t) 

Note that, recent studies (see for example [6]) show that 
(cr — r(f)) is around 6.5 days. 

III. Data fitting 

In this section we provide some numerical experiments 
based on data from Guinea, Sierra-Leon and liberia. We 
consider the cumulative number of infectious cases and 
deaths denoted by C(t) and respectively. They can 

be calculated as 

t 

C(t + 1) =^2x(s - d); (5) 

s—0 

t n 

D(t + 1) = EE auj p (i) x(s — d — i). (6) 

S—0 2—0 

Here Ooj p (i) is the death rate of infectious population in 
generation x(s d — i) and n is a large number. Parameters 
of the gamma distribution function cu(i) are taken from [15] 
where 

lj p (x) = a = 10, b= 1.3333 (7) 

I [a) 

with mean value 7.5. Moreover, we set d = 7 and n = 35 
as in [7]. 

In the considered model 0 there are only three 
parameters a and fj (constants) and a piece-wise constant 
control function r(f) that need to be optimized to fit data - 
the total number of infectious population and deaths. The 
aim here is to show that there exists a control r(f) such that 
the corresponding trajectory x(t) fits data well. 



TABLE I 

Results of best fits: the (effective) reproduction 

NUMBERS Rk AND AVERAGE TIMES TO HOSPITALIZATION Tk (IN 
DAYS) FOR DIFFERENT INTERVALS Afc, k = 1, 2, 3, 4. THE 
OPTIMAL VALUES FOR a AND /? ARE ALSO PROVIDED; THEY 
ARE CONSTANT FOR A WHOLE PERIOD 


Country 

a 

P 

Rl (Tl) 

R -2 (T2) 

R3 O3) 

i?4 (T 4 ) 

Gui. 

0.66 

0.265 

0.79 (3) 

1.31 (5) 

1.06 (4) 

0.79 (3) 

S.-L. 

0.32 

0.274 

1.36 (5) 

1.36 (5) 

1.09 (4) 

0.82 (3) 

Lib. 

0.46 

0.294 

1.17 (4) 

1.46 (5) 

0.88 (3) 

0.88 (3) 


Guinea 



We consider three consequent intervals A*. = [T k , T k+1 ] 
(k = 1, • - -, 4) for each country and find optimal values 
a, /3 and Tk ( k = 1, • • ■,4) where r(f) = Tk, Vi £ A^. 
The results are presented in Table [I] The last time point T 5 
is 01-Mar-2015. The values of T \ , 7j, 7 3, 7 4 are as follows: 
22-March, 23-May, 20-July and 04-Dec-2014 for Guinea; 27- 
May, 20-June, 20-August and 04-Dec-2014 for Sierra Leone; 
and 16-June, 20-July, 07-Sept and 04-Dec-2014 for Liberia. 
Each interval A/, has its own reproduction number Ilk that 
defines the shape of the best fits presented in Figure [I] 

Data were retrieved from the WHO website [17] for the 
cumulative numbers of clinical cases (confirmed, probable 
and suspected) collected till 1 March 2015. The global 
optimization algorithm DSO in Global And Non-Smooth 
Optimization (GANSO) library [18], [19] is applied for 
solving optimization problems in this section as well as in 
Section [V| 

The results obtained show that the estimated values of a 
and /3 can be used for future time intervals while considering 
rasa dynamic parameter that defines the spread of infection. 
This naturally leads to optimal control problems that are 
considered in the next section. 

IV. OPTIMAL DISTRIBUTION OF BED 
CAPACITIES 

Denote by B(t) the number of beds (capacity of hospitals) 
at time t. It is an increasing and piece-wise constant function 
where jumps are at the points T- 1 , % = 1, 2. ■■■,(/. We refer 
to [ 6 ] (Fig 1) for an example of B(t) in Liberia between 
June-September 2014. Note that q > p\ that is the number 
of points r lf is larger than T r 


Sierra-Leon 



Liberia 



Fig. 1. The best fits for the cumulative numbers of infected cases and 
deaths in Guinea, Sierra Leone and Liberia by considering three parameters 
a, ft and 4 subintervals with different values r*., k = 1,2, 3, 4 (for the 
values see Table |TJ. The lines represent the best fits, red and black circles 
represent the data 

down. Taking into account this observation, we use the 
hospitalization rate in the definition of the objective function 
(to be minimized) given by 


We assume that the initial period of epidemic is long 
enough to estimate parameters a, /3 and to provide some 
future scenarios depending on r(t). 


T 

x(T)+K-J2 

t =1 


h(t : r, x) 
B(t) 


where K is constant. 


A. Objective function 

Given process (r, x) we define hospitalization rate as 

HospRatejt) = ■ ( 8 ) 

B(t) 

Note that the number of infectious cases (especially during 
high growth) may exceed the number of available beds, 
while it becomes quite low when the infection is slowing 


B. Constraints in terms of costs 

We consider the following three functions that will be used 
to formulate cost constraints. 

• cb(A&) - costs associated to building A b additional 
beds; 

• cs{h ) - costs required for servicing h patients in hos¬ 
pitals; 

• Ci[Ah) - costs required for Ah infectious cases before 
hospitalization. 




































Assuming that F(t) is the total available funds, we can 
formulate cost constraints as 

c s (A6(t)) +c s {h(t)) +c/(A h(t)) < F(t), Vt. 

Here A6(f) = b(t + 1) — 6(f) and A h(t) = h(t + 1 : 
r, x) — h(t : r,x). 

C. Multi regions 

Now consider m regions/countries and assume there is 
an inter-transmission of infections between them. Denote 
by x r (t) the number of infected cases in country r and 
assume that the transmission of infection from country i 
to r is given by the coefficient f3 ir ([13]). In this case, by 
considering corresponding average times of isolation r r (f) 
in each country r, we have the following system 

x r {t + 1) = /3ir Eliott 1 ~ aiu(i))xi(t -d-i ) 

+/?2 r ES^E 1 ~ “2 U>(i)) X 2 (t -d-i) 

+/3mr El=O t)_1 ( 1 ~ a mW(t)) X m (t - d - l) 

where /3j r is the infections generated from country i. One 
would expect fi rr to be much larger than (3 ir for i ^ r 
and in a special case /3i r = 0. The number of hospitalized 
population at time t is defined by 0; that is 


h r (t: r r ,x r ) = (l—a r ut(i))x(t—d—i), r = 

»=Tr(t) 

We assume that the death rates a r and the coefficients B vr 
are estimated from the initial data and they are constant over 
the whole period of the epidemic. Then, we can consider 
the problem of optimal distribution of available new beds 
between regions formulated below. 

Problem 1 (Optimal distribution of bed capacities): 

Given initial data (a r , fdir and x r (t),t < 1), total bed 
capacity B(t) and future scenario for control functions 
T r (t),r = I,- — ,to, find increasing piece-wise constant 
functions 6 r (t), r = 1 , • • •, to , for the problem 


Minimize : 


E 


x r (T) + K 


T 

E 

t=i 


h r (t ■ TV, x r ) 

b r (t) 


subject to : &i(t) H— • + b m (t) < B(t ), Vi. 


In the next problem we take into account the cost 
constraints: 

Problem 2 (Optimal distribution of bed capacities 
under cost constraints): Given initial data (a r , B lr and 
x r (t),t < 1), total budget function F(t) and future scenario 
for control functions Ty(f),r = 1 ,---,to, find increasing 
piece-wise constant functions b r (t),r = 1 for the 

problem 


Minimize : 


-(- T) + I<-Y.] 


h r (t . T r , x r 

b r (t) 


subject to : 6 i(i) H— • + 6 m (f) < B(t), Vi; 


c B (AS(i)) + ^2 [ c s(K{t)) + c/(A h r (t))} < F(t), Vi. 

r—1 

Feasible processes. Given initial data (a r> B, r and 
x r {t),t < 1) and B(t)) consider the trajectory x = 

(xi,---,x m ) corresponding to r = (ri,---,r m ). Denote 
also b = (bi, • • •, b m ), where each b r stands for the bed 
capacity function b r (t) in region r. 

We call process (r, b, x) feasible if all the constraints of 
the problem under consideration hold and the hospitalization 
rates are less than 1 ; that is, 


h r (t : r r , x r ) < b r (t), Vi, r = 1, • • •, to. 

The meaning of feasible processes can be explained as 
follows. If (r, b, x) is feasible then the number of required 
beds and the resources needed for isolation are sufficient at 
every time point f in order to keep the average times of 
isolation at level r = (ti, • • •, r m ). Thus, the corresponding 
effective reproduction numbers r can be considered as upper 
bounds (the actual effective reproduction numbers might be 
even lower). Therefore, a feasible process determines in 
some sense the best use of given resources to achieve the 
“guaranteed lowest” number of infectious cases. 


V. NUMERICAL EXPERIMENTS 

In this section we provide an example on a synthetic 
data set to demonstrate how the problems formulated above 
can be used for controlling the spread of infection. In this 
example there are two regions (to = 2 ) and for the sake 
of simplicity we assume that there is no transmission of 
infection between these regions (i.e. = /3 2,1 = 0 ). 

Moreover, we consider only Problem 1; that is, costs related 
to bed building (cb), services ( cs ) and before isolation (c/) 
assumed to be sufficient in all cases. 


Initial data. We assume the time interval is [1, T\ = 

[1.150] ; where [1,100] is an initial (past) period and 

[101.150] is the future/planning period for our optimal 
control problem. 

The number of initial (i.e. t < 0) infected cases is 2 in 
both regions. The set of possible values for r is {3,4,5} 
(as in the case of data fitting in Section m- We generate 
synthetic data - Xi(t), x 2 (t) for t £ [1,100] by setting 

• Region 1: a = 0.6, /3ii = 0-30 and T±(t) = 4, Vt £ 
[1,50], Ti(f) = 5,Vt £ [51,100]; 

• Region 2: a = 0.6, B 22 = 0.28 and 72(f) = 4,Vt £ 
[1,50], 72 (f) = 5,Vt £ [51,100], 

As in [7], oj as a gamma distribution function with 
mean value 7.5 defined by 0. According to formula 0 
corresponding effective reproduction numbers are 










• Region 1: R = 0.90,1.20 and 1.48 for r = 3,4 and 5, 
respectively; 

• Region 2: R = 0.84,1.12 and 1.38 for r = 3,4 and 5, 
respectively. 

Thus, the initial functions Xi(t),X 2 (t) for t £ [1,100] have 
effective reproduction numbers R = 1.20 and 1.48 for x\(t) 
on [1, 50] and [51,100], respectively; R = 1.12 and 1.38 for 
X 2 (f) on [1,50] and [51,100], respectively. 

We will consider Problem 1 on the interval [101,150]. 
Both controls Ti(t) will be assumed to be constant: iy(t) = 
Ti , Vt £ [101,150], i = 1,2. Values ti and r 2 will be used 
to describe future possible scenarios. 

The initial number of beds are 6i(100) = 126 and 
62 ( 100 ) = 60. We assume that new beds will be created 
at times T-f = 101, T| = 108, T| = 115 and T| = 122; 
corresponding numbers of additional beds will be denoted 
by Abi, i = 1,2,3,4. Optimal control problem aims to 
distribute these additional beds between the regions. 

We introduce a new variable - A i that denotes part of A 6 ,; 
considered for the first region, the remaining part (1 —Ai)A 6 i 
for the second region. 

Under these assumptions, Problem 1 can formulated as 
follows: 

Minimize (AljA2;A3jA4) Zi(150) + £ 2 (150) 
hi(t : t 1 ,x 1 ) h 2 {t : t 2 ,x 2 ) ~ 

61(f) 62(f) 

subject to : A k £ [0,1], k = 1,2, 3 , 4 ; 

x r (t + 1) = P rr (1 — aco(i)) x r (t — d — i), r = 1, 2; 

i =0 

|{Tj><t; j —1,2,3,4}| 

61(f) = 6i(100) + Y A*-A6i; 

i=i 

\{Tj<t\ j= 1,2,3,4}| 

6 2 (f) = 62(100)+ Y ( 1 -A 0 -A 6 * 

i=1 

Here K = 100 and h r (t : r r ,x r ) is calculated by 0 by 
setting cr = 6 and d = 6 . 

We can consider different scenarios depending on T \, r 2 
and additional beds A 6 &, k = 1,2,3,4. Below we will 
provide three scenarios where we set 

A 6 i = 350, A 6 2 = 300, A 6 3 = 100 , A 6 4 = 20 

and change the values of r 4 , T 2 as 3, 4 and 5 (that is, all 
possible values used in the data fitting problem in Section 
0 - The aim here is to compare corresponding optimal 
distributions of bed capacities. We recall that additional beds 
are introduced (weekly) at T-f = 101, T| = 108, T 3 b = 115 
and T| = 122. The cumulative number of infected cases at 


the start f = 100 are: 1259 in Region 1 and 675 in Region 2. 

Case 1: t± = T 2 = 3. 

The summary of the optimal solution obtained is provided 
below. 

• The optimal distribution of additional beds: 

Region 1 : 192.3 183.4 72.9 20 

RegionZ : 157.7 116.6 27.1 0 

Total : 350 300 100 20 

• Cumulative number of infected cases at the end (f = 
150) of planning period: 2951 in Region 1 and 1289 in 
Region 2. 

• The average number of bed occupancy over the time 
interval [100,150] is 0.45 for Region 1 and 0.29 for 
Region 2. 

• The maximum occupancy rates are: 0.83 (that is, av¬ 
erage 0.83 patient per bed) in Region 1 and 0.55 in 
Region 2; that is, the demand for hospital beds is met 
at every point f £ [100,150]. 

Thus the solution obtained is feasible. 

Case 2: r 4 = r 2 = 4. 

Optimal solution obtained and some relevant parameters 
are: 

• The optimal distribution of additional beds: 

Region! : 192.0 183.4 74.5 20 

Region2 : 158.0 116.6 25.5 0 

Total : 350 300 TOO W 

• Cumulative number of infected cases at the end (f = 
150) of planning period: 5846 in Region 1 and 2259 in 
Region 2. 

• The average number of bed occupancy over the time 
interval [100,150] is 0.65 for Region 1 and 0.40 for 
Region 2. 

• The maximum occupancy rates are: 0.96 (that is, av¬ 
erage 0.96 patient per bed) in Region 1 and 0.52 in 
Region 2. 

Again, the demand for hospital beds is met at every point 
f £ [100,150] and accordingly the optimal solution obtained 
is feasible. 

Case 3: r 4 = r 2 = 5. 

Optimal solution obtained and some relevant parameters 
are: 

• The optimal distribution of additional beds: 

Regionl : 191.6 183.3 75 20 

Region2 : 158.4 116.7 25 0 

Total : 350 300 100 20 

• Cumulative number of infected cases at the end it = 
150) of planning period: 11585 (Region 1) and 4163 
(Region 2). 


150 

+K- Y, 









TABLE II 

The optimal distribution of additional beds 


n 

T2 


3 

4 

Region 1 

204.7 

183.9 

9.1 

0 



Region 2 

145.3 

116.1 

90.9 

20 

4 

5 

Region 1 

206.2 

188.2 

31.6 

0 



Region 2 

143.8 

111.8 

68.4 

20 

4 

3 

Region 1 

179.3 

223.3 

100 

20 



Region 2 

170.7 

76.7 

0 

0 

5 

4 

Region 1 

177.1 

208.5 

100 

20 



Region 2 

172.9 

91.5 

0 

0 


• The average number of bed occupancy over the time 
interval [100,150] is 0.83 for Region 1 and 0.51 for 
Region 2. 

• The maximum occupancy rates are: 1.91 (that is, 1.91 
patient per bed) in Region 1 and 1.05 in Region 2. 

Therefore this solution is not feasible as the bed capacities 
are not enough for isolation all infected individuals. 

Comparting the results in Cases 1-3, where t\ = T 2 , we 
observe that, the optimal distributions of bed capacities are 
almost the same although the solution obtained in Case 3 
is even not feasible. In table [II] we also provide the results 
obtained by assuming that the rate of increase in one region 
is greater than the other one. The results for tt = T 2 + I (that 
is, Ti = 3, T 2 = 4 and ti = 4, T 2 = 5) display quite similar 
optimal bed distributions. We observe the same situation for 

r 2 = Ti + 1. 

Summarizing these results we note that the optimal con¬ 
trol problem considered can provide quite “robust” optimal 
distributions of new bed capacities across the regions under 
each of the assumptions t\ = T 2 , tt > T 2 and n < 72 . 
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